home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SPARSE.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  83 lines
  1. PROCEDURE sparse(b: glnarray; n: integer; VAR x: glnarray; VAR rsq: real);
  2. (* Programs using routine SPARSE must define the type
  3. TYPE
  4.    glnarray = ARRAY [1..n] OF real;
  5. in the main routine. They must also provide two routines,
  6. PROCEDURE asub(x: glnarray; VAR y: glnarray; n: integer);
  7. and
  8. PROCEDURE atsub(x: glnarray; VAR z: glnarray; n: integer);
  9. which calculate A*x and (A transpose)*x *)
  10. LABEL 1,99;
  11. CONST
  12.    eps=1.0e-6;
  13. VAR
  14.    j,iter,irst: integer;
  15.    rp,gg,gam,eps2,dgg,bsq,anum,aden: real;
  16.    g,h,xi,xj: glnarray;
  17. BEGIN
  18.    eps2 := n*sqr(eps);
  19.    irst := 0;
  20. 1:   irst := irst+1;
  21.    asub(x,xi,n);
  22.    rp := 0.0;
  23.    bsq := 0.0;
  24.    FOR j := 1 TO n DO BEGIN
  25.       bsq := bsq+sqr(b[j]);
  26.       xi[j] := xi[j]-b[j];
  27.       rp := rp+sqr(xi[j])
  28.    END;
  29.    atsub(xi,g,n);
  30.    FOR j := 1 TO n DO BEGIN
  31.       g[j] := -g[j];
  32.       h[j] := g[j]
  33.    END;
  34.    FOR iter := 1 TO 10*n DO BEGIN
  35.       asub(h,xi,n);
  36.       anum := 0.0;
  37.       aden := 0.0;
  38.       FOR j := 1 TO n DO BEGIN
  39.          anum := anum+g[j]*h[j];
  40.          aden := aden+sqr(xi[j])
  41.       END;
  42.       IF (aden = 0.0) THEN BEGIN
  43.          writeln('pause in routine SPARSE');
  44.          writeln('very singular matrix'); readln
  45.       END;
  46.       anum := anum/aden;
  47.       FOR j := 1 TO n DO BEGIN
  48.          xi[j] := x[j];
  49.          x[j] := x[j]+anum*h[j]
  50.       END;
  51.       asub(x,xj,n);
  52.       rsq := 0.0;
  53.       FOR j := 1 TO n DO BEGIN
  54.          xj[j] := xj[j]-b[j];
  55.          rsq := rsq+sqr(xj[j])
  56.       END;
  57.       IF ((rsq = rp) OR (rsq <= bsq*eps2)) THEN GOTO 99;
  58.       IF (rsq > rp) THEN BEGIN
  59.          FOR j := 1 TO n DO BEGIN
  60.             x[j] := xi[j]
  61.          END;
  62.          IF (irst >= 3) THEN GOTO 99;
  63.          GOTO 1
  64.       END;
  65.       rp := rsq;
  66.       atsub(xj,xi,n);
  67.       gg := 0.0;
  68.       dgg := 0.0;
  69.       FOR j := 1 TO n DO BEGIN
  70.          gg := gg+sqr(g[j]);
  71.          dgg := dgg+(xi[j]+g[j])*xi[j]
  72.       END;
  73.       IF (gg = 0.0) THEN GOTO 99;
  74.       gam := dgg/gg;
  75.       FOR j := 1 TO n DO BEGIN
  76.          g[j] := -xi[j];
  77.          h[j] := g[j]+gam*h[j]
  78.       END
  79.    END;
  80.    writeln('pause in routine SPARSE');
  81.    writeln('too many iterations'); readln;
  82. 99:   END;
  83.